#ifndef _imageTypes_h_
#define _imageTypes_h_

// imageTypes.h
//
// Definition of our supported image types

#ifdef WIN32
#pragma warning(disable:4786) // Disable warning regarding long symbol names caused by templates
#pragma warning(disable:4244) // Don't show possible conversion error for op+= and related
#endif

#include "geometry.h"
#include <iostream>

namespace simple{
    //
    //LimitInfo
    //
    // A generic image processing function needs to know certain
    // information about data types. For example, the minimum value
    // for a type is something which is difficult to decide at runtime.
    // A TypeInfo object should be defined for each datatype. This
    // allows proper rounding, clamping, etc. to occur.
    //
    // One object is constructed per type. This object is similar to the
    // numeric_limits<> defined in <limits>, except we only care about 
    // min() and max().

    template<class T> class LimitInfo
    {
    public:
        static LimitInfo<T> sType;

        T minValue;
        T maxValue;

    private:
        LimitInfo (T min, T max) : minValue (min), maxValue (max) {}
    };


    // Perform limit checking so pixel values don't wrap
    // Usage: D dst = Limit<D> (src); 
    // Use this instead of 'dst= src'
    // Some specializations are found elsewhere in this file
    template<class D, class S> D Limit (const S& src)
    {
        // 边界内取原值，边界外取边界值。
        return src < LimitInfo<D>::sType.minValue ? LimitInfo<D>::sType.minValue : 
            (src > LimitInfo<D>::sType.maxValue ? LimitInfo<D>::sType.maxValue : static_cast<D>(src));
    }



    //
    // Generic Functions
    //
    // The image processing functions frequently call functions
    // like this to perform basic operations. These do so with
    // clamping.
    // (clamping : 时期可以象其它数据类型一样使用，
    // 并消除由于溢出问题而近期的不必要的像素包装行为。)
    //
    // We would very much like to define a generic function template
    // that looks like,
    //
    //    template<class D, class S1, class S2>
    //    D add2 (const S1& s1, const S2& s2)
    //    {
    //      return static_cast<D>(s1) + s2;
    //    }
    //
    // However, we aren't so fortunate, because we also want to 
    // define a function that looks like this,
    //
    //    template<class D, class S1, class S2>
    //    D add2 (const ClampedTmpl<S1>& s1, const ClampedTmpl<S2>& s2)
    //    {
    //      return Limit<D> (s1.val + s2.val);   // d = s1 + s2;
    //    }
    //
    // The problem is that these create an ambiguity when a template
    // of the latter type (ie. ClampedTmpl<>) is instantiated. See
    // the C++ standard, "Partial Ordering of Function Templates" 
    // (section 5.5.2) and "Template Argument Deduction" (section 8.2).
    //
    // The behavior on various compilers is erratic. To avoid this problem,
    // we explicitly define the common forms we'll need and let the compiler
    // generate the other ones as needed. This avoid the ambiguity error.

    #if 0
    // early gcc can handle this without problem. For now, we'll define
    // these by hand.

    template<class D, class S1, class S2>
    D add2 (const S1& s1, const S2& s2)
    {
        return static_cast<D>(s1) + s2;
    }


    template<class D, class S1, class S2>
    D sub2 (const S1& s1, const S2& s2)
    {
        return static_cast<D>(s1) - s2;
    }

    template<class D, class S1, class S2>
    D mul2 (const S1& s1, const S2& s2)
    {
        return static_cast<D>(s1) * s2;
    }

    template<class D, class S1, class S2>
    D div2 (const S1& s1, const S2& s2)
    {
        return static_cast<D>(s1) / s2;
    }

    template<class D, class S>
    D scale (S& s1, float scaling)
    {
        return static_cast<D> (scaling * static_cast<float>(s1));
    }


    #else
    // We define manual versions. If you encounter errors with missing versions
    // of add2<>, sub2<>, mul2<>, div2<>, or scale<>, chances are you will need
    // to add more entries here.

    template<class D> D add2 (Pel32s s1, Pel32s s2)
    { return static_cast<D>(s1 + s2);}

    template<class D> D sub2 (Pel32s s1, Pel32s s2)
    { return static_cast<D>(s1 - s2);}

    template<class D> D mul2 (Pel32s s1, Pel32s s2)
    { return static_cast<D>(s1 * s2);}

    template<class D> D div2 (Pel32s s1, Pel32s s2)
    { return static_cast<D>(s1 / s2);}

    template<class D> D scale (Pel32s s1, float scaling)
    { return static_cast<D> (scaling * s1);}

    template<class D> D scale (Pel8 s1, float scaling)
    { return static_cast<D> (scaling * s1);}

    template<class D> D scale (Pel16 s1, float scaling)
    { return static_cast<D> (scaling * s1);}

    template<class D> D scale (double s1, float scaling)
    { return static_cast<D> (scaling * s1);}

    #endif

    //
    // Clamped pixel type 
    //
    // This list of supported operators is not necessarily 
    // complete, but it is sufficient for our class right now.
    //
    // Everything is defined in the header file. The definitions are simple
    // and win32 (at least VC6) has problems when you split the declaration
    // from the definition (and use the 'export' keyword).
    template<class T> class ClampedTmpl
    {
    public:
        T val;

        operator T () const { return val;}

        // constructors
        ClampedTmpl ()    : val(0) {}

        ClampedTmpl (T v) : val(v) {}

        ClampedTmpl (const ClampedTmpl& src)
        {
            val = Limit<T> (src.val);
        }

        template<class T1>
        ClampedTmpl (const ClampedTmpl<T1>& src)
        {
            val = Limit<T> (src.val);
        }

        ClampedTmpl<T>& operator= (const ClampedTmpl& src)
        {
            val = Limit<T> (src.val);
            return *this;
        }

        /*
        template<class T1>
        ClampedTmpl<T>& operator= (const ClampedTmpl<T1>& src)
        {
            val = Limit<T> (src.val);
            return *this;
        }
        */

        /*
        ClampedTmpl& operator= (const T& c)
        {
            val = c;
            return *this;
        }
        */

        ClampedTmpl<T>& operator+= (const ClampedTmpl<T>& s)
        {
            val = Limit<T> (val + s.val);
            return *this;
        }

        ClampedTmpl<T>& operator-= (const ClampedTmpl<T>& s)
        {
            val = Limit<T> (val - s.val);
            return *this;
        }

        ClampedTmpl<T>& operator*= (const ClampedTmpl<T>& s)
        {
            val = Limit<T> (val * s.val);
            return *this;
        }

        ClampedTmpl<T>& operator/= (const ClampedTmpl<T>& s)
        {
            val = Limit<T> (val / s.val);
            return *this;
        }

        template <class T2>
        ClampedTmpl<T> operator+ (const ClampedTmpl<T2>& s2)
        {
            ClampedTmpl<T> dst;
            dst.val = Limit<T> (val + s2.val);
            return dst;
        }

    };


    //
    // typedef
    //
    // Common types our applications might use

    typedef ClampedTmpl<Pel8>   ClampedPel8;
    typedef ClampedTmpl<Pel16>  ClampedPel16;
    typedef ClampedTmpl<Pel32>  ClampedPel32;
    typedef ClampedTmpl<Pel32s> ClampedPel32s;   // This is signed!!



    //
    // ClampedTmpl<T> non-member
    //
    // bool operator== (const ClampedTmpl<T>& s1, const ClampedTmpl<T>& s2);
    // bool operator== (const ClampedTmpl<T>& s1, const T& s2);
    // bool operator!= (const ClampedTmpl<T>& s1, const ClampedTmpl<T>& s2);
    // bool operator!= (const ClampedTmpl<T>& s1, const T& s2);
    //
    // ClampedTmpl<T1> operator+ (const ClampedTmpl<T1>& s1, const T2& s2);
    // ClampedTmpl<T1> operator- (const ClampedTmpl<T1>& s1, const T2& s2);
    // ClampedTmpl<T1> operator* (const ClampedTmpl<T1>& s1, const T2& s2);
    // ClampedTmpl<T1> operator/ (const ClampedTmpl<T1>& s1, const T2& s2);
    //
    // std::ostream& operator<< (std::ostream& stream, const ClampedTmpl<T>& val);


    template<class T>
    bool operator== (const ClampedTmpl<T>& s1, const ClampedTmpl<T>& s2)
    {
        return s1.val == s2.val;
    }

    template<class T>
    bool operator!= (const ClampedTmpl<T>& s1, const ClampedTmpl<T>& s2)
    {
        return s1.val != s2.val;
    }


    template<class T>
    bool operator== (const ClampedTmpl<T>& s1, const T& s2)
    {
       return s1.val == s2;
    }

    template<class T>
    bool operator!= (const ClampedTmpl<T>& s1, const T& s2)
    {
        return s1.val != s2;
    }


    template<class T1, class T2>
    ClampedTmpl<T1> operator+ (const ClampedTmpl<T1>& s1, const T2& s2)
    {
        ClampedTmpl<T1> dst;
        dst.val = Limit<T1> (s1.val + s2);
        return dst;
    }

    template<class T1, class T2>
    ClampedTmpl<T1> operator- (const ClampedTmpl<T1>& s1, const T2& s2)
    {
        ClampedTmpl<T1> dst;
        dst.val = Limit<T1> (s1.val - s2);
        return dst;
    }

    template<class T1, class T2>
    ClampedTmpl<T1> operator* (const ClampedTmpl<T1>& s1, const T2& s2)
    {
        ClampedTmpl<T1> dst;
        dst.val = Limit<T1> (s1.val * s2);
        return dst;
    }

    template<class T1, class T2>
    ClampedTmpl<T1> operator/ (const ClampedTmpl<T1>& s1, const T2& s2)
    {
        ClampedTmpl<T1> dst;
        dst.val = Limit<T1> (s1.val / s2);
        return dst;
    }


    template<class T>
    std::ostream& operator<< (std::ostream& stream, const ClampedTmpl<T>& val)
    {
        stream << val.val;
        return stream;
    }

    // We want to display the pixel value as a number, not as a character
    // Even though this is an unsigned quantity, some compilers make mistakes
    template<> std::ostream& operator<< (std::ostream& stream, const ClampedTmpl<Pel8>& rgb);




    //
    // ClampedTmpl<T> mixed
    //
    // add2<D> (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // sub2<D> (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // mul2<D> (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // div2<D> (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // scale<D> (const ClampedTmpl<T1>& s1, float scaling);
    //
    // ClampedTmpl<T1> operator+ (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // ClampedTmpl<T1> operator- (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // ClampedTmpl<T1> operator* (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // ClampedTmpl<T1> operator/ (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    //
    // ClampedTmpl<T1>& operator+= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // ClampedTmpl<T1>& operator-= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // ClampedTmpl<T1>& operator*= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);
    // ClampedTmpl<T1>& operator/= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2);


    // Limit<> specializations
    template<> ClampedPel8 Limit<ClampedPel8, ClampedPel32s> (const ClampedPel32s& src);

    template<> Pel8   Limit<Pel8,   Pel32s> (const Pel32s& src);
    template<> Pel8   Limit<Pel8,   Pel8>   (const Pel8& src);
    template<> Pel32s Limit<Pel32s, Pel32s> (const Pel32s& src);
    template<> Pel32s Limit<Pel32s, Pel32>  (const Pel32& src);

    template<> ClampedPel8 Limit<ClampedPel8, ClampedPel8> (const ClampedPel8& src);

    template<> ClampedPel32s Limit<ClampedPel32s, ClampedPel32s> (const ClampedPel32s& src);


    template<class D, class S1, class S2>
    D add2 (const ClampedTmpl<S1>& s1, const ClampedTmpl<S2>& s2)
    {
        return Limit<D> (s1.val + s2.val);   // d = s1 + s2;
    }

    template<class D, class S1, class S2>
    D sub2 (const ClampedTmpl<S1>& s1, const ClampedTmpl<S2>& s2)
    {
       return Limit<D> (s1.val - s2.val);   // d = s1 - s2;
    }

    template<class D, class S1, class S2>
    D mul2 (const ClampedTmpl<S1>& s1, const ClampedTmpl<S2>& s2)
    {
       return Limit<D> (s1.val * s2.val);   // d = s1 * s2;
    }

    template<class D, class S1, class S2>
    D div2 (const ClampedTmpl<S1>& s1, const ClampedTmpl<S2>& s2)
    {
       return Limit<D> (s1.val / s2.val);   // d = s1 / s2;
    }

    template<class D, class S1>
    D scale (const ClampedTmpl<S1>& s1, float scaling)
    {
       return Limit<D> (scaling * static_cast<float>(s1.val));   // d = k1 * s1
    }


    template<class T1, class T2>
    ClampedTmpl<T1> operator+ (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        ClampedTmpl<T1> dst = add2<ClampedTmpl<T1>, ClampedTmpl<T1>, ClampedTmpl<T2> >(s1, s2);
        return dst;
    }

    template<class T1, class T2>
    ClampedTmpl<T1> operator- (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        ClampedTmpl<T1> dst = sub2<ClampedTmpl<T1>, ClampedTmpl<T1>, ClampedTmpl<T2> >(s1, s2);
        return dst;
    }

    template<class T1, class T2>
    ClampedTmpl<T1> operator* (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        ClampedTmpl<T1> dst = mul2<ClampedTmpl<T1>, ClampedTmpl<T1>, ClampedTmpl<T2> >(s1, s2);
        return dst;
    }

    template<class T1, class T2>
    ClampedTmpl<T1> operator/ (const ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        ClampedTmpl<T1> dst = div2<ClampedTmpl<T1>, ClampedTmpl<T1>, ClampedTmpl<T2> >(s1, s2);
        return dst;
    }


    template<class T1, class T2>
    ClampedTmpl<T1>& operator+= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        s1.val = Limit<T1> (s1.val + s2.val);   // d = s1 + s2;
        return s1;
    }

    template<class T1, class T2>
    ClampedTmpl<T1>& operator-= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        s1.val = Limit<T1> (s1.val - s2.val);   // d = s1 - s2;
        return s1;
    }

    template<class T1, class T2>
    ClampedTmpl<T1>& operator*= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        s1.val = Limit<T1> (s1.val * s2.val);   // d = s1 * s2;
        return s1;
    }

    template<class T1, class T2>
    ClampedTmpl<T1>& operator/= (ClampedTmpl<T1>& s1, const ClampedTmpl<T2>& s2)
    {
        s1.val = Limit<T1> (s1.val / s2.val);   // d = s1 / s2;
        return s1;
    }



    //
    // RGBTmpl<T>
    //
    // A generic rgb triplet. I originally wanted to name it RGB but
    // this runs into problems with macros defined on some platforms
    // (like win32). The data members are left public so it acts like
    // a struct. There are no virtual functions so the size
    // of this object matches how we want memory to be written in
    // memory. Most operators are written as non-member functions.

    template<class T> class RGBTmpl
    {
    public:
        T red;
        T green;
        T blue;

        RGBTmpl ()              : red(0), green(0), blue(0) {}
        explicit RGBTmpl (T v)  : red(v), green(v), blue(v) {}
        RGBTmpl (T r, T g, T b) : red(r), green(g), blue(b) {} 

        T ts1 () const { return red;}
        T ts2 () const { return green;}
        T ts3 () const { return blue;}
        // Treat as generic tristimulas values

        RGBTmpl (const RGBTmpl& s)
        {
            red   = s.red;
            green = s.green;
            blue  = s.blue;
        }
        // Copy constructor

        RGBTmpl<T>& operator= (const RGBTmpl& src)
        {
            red   = src.red;
            green = src.green;
            blue  = src.blue;
            return *this;
        }

        template<class T1> RGBTmpl (const RGBTmpl<T1>& s)
        {
            red   = Limit<T> (s.red);
            green = Limit<T> (s.green);
            blue  = Limit<T> (s.blue);
        }
        // Copy constructor

        template<class T1> RGBTmpl<T>& operator= (const RGBTmpl<T1>& src)
        {
            red   = Limit<T> (src.red);
            green = Limit<T> (src.green);
            blue  = Limit<T> (src.blue);
            return *this;
        }

        RGBTmpl& operator= (const T& c)
        {
            red   = c;
            green = c;
            blue  = c;
            return *this;
        }

        operator T () const
        { return (red + green + blue) / 3;}
        // Conversion to monochrome (averaging)

        // If you want the NTSC definition, use this
        //  double v = .299 * red + .587 * green + .114 * blue;
        //  return static_cast<T> (v);


        // A conversion operator is not needed because it does 
        // the same job as this:
        //  template<class T1> RGBTmpl (const RGBTmpl<T1>& s)
        //
        // Most compilers we tested did not complain about it.
        //
        //   template<class T2> operator RGBTmpl<T2> () const
        //   {
        //      RGBTmpl<T2> dst;
        //      dst.red   = Limit<T> (red);
        //      dst.green = Limit<T> (green);
        //      dst.blue  = Limit<T> (blue);
        //      return dst;
        //    }


        RGBTmpl<T>& operator+= (const RGBTmpl<T>& s)
        {
            red   = Limit<T> (red   + s.red);
            green = Limit<T> (green + s.green);
            blue  = Limit<T> (blue  + s.blue);
            return *this;
        }

        RGBTmpl<T>& operator-= (const RGBTmpl<T>& s)
        {
            red   = Limit<T> (red   - s.red);
            green = Limit<T> (green - s.green);
            blue  = Limit<T> (blue  - s.blue);
            return *this;
        }

        RGBTmpl<T>& operator*= (const RGBTmpl<T>& s)
        {
            red   = Limit<T> (red   * s.red);
            green = Limit<T> (green * s.green);
            blue  = Limit<T> (blue  * s.blue);
            return *this;
        }

        RGBTmpl<T>& operator/= (const RGBTmpl<T>& s)
        {
            red   = Limit<T> (red   / s.red);
            green = Limit<T> (green / s.green);
            blue  = Limit<T> (blue  / s.blue);
            return *this;
        }


        RGBTmpl<T>& operator+= (const T& s)
        {
            red   = Limit<T> (red   + s);
            green = Limit<T> (green + s);
            blue  = Limit<T> (blue  + s);
            return *this;
        }

        RGBTmpl<T>& operator-= (const T& s)
        {
            red   = Limit<T> (red   - s);
            green = Limit<T> (green - s);
            blue  = Limit<T> (blue  - s);
            return *this;
        }

        RGBTmpl<T>& operator*= (const T& s)
        {
            red   = Limit<T> (red   * s);
            green = Limit<T> (green * s);
            blue  = Limit<T> (blue  * s);
            return *this;
        }

        RGBTmpl<T>& operator/= (const T& s)
        {
            red   = Limit<T> (red   / s);
            green = Limit<T> (green / s);
            blue  = Limit<T> (blue  / s);
            return *this;
        }
    };




    //
    // typedef
    //
    typedef RGBTmpl<Pel8>   RGB;
    typedef RGBTmpl<Pel32>  RGBPel32;
    typedef RGBTmpl<Pel32s> RGBPel32s;   // This is signed!!


    //
    // RGBTmpl<T> non-member
    //
    //   operator== (const RGBTmpl<T>& s1, const RGBTmpl<T>& s2)
    //   operator!= (const RGBTmpl<T>& s1, const RGBTmpl<T>& s2)
    //
    //   operator+ (const RGBTmpl<T1>& s1, const T2& s2)
    //   operator- (const RGBTmpl<T1>& s1, const T2& s2)
    //   operator* (const RGBTmpl<T1>& s1, const T2& s2)
    //   operator/ (const RGBTmpl<T1>& s1, const T2& s2)
    //
    //   operator+ (const RGBTmpl<T1>& s1, const double& s2)
    //   operator- (const RGBTmpl<T1>& s1, const double& s2)
    //   operator* (const RGBTmpl<T1>& s1, const double& s2)
    //   operator/ (const RGBTmpl<T1>& s1, const double& s2)

    //   operator<< (std::ostream& stream, const RGBTmpl<T>& rgb)

    template<class T>
    bool operator== (const RGBTmpl<T>& s1, const RGBTmpl<T>& s2)
    {
        return (s1.red == s2.red) && (s1.green == s2.green) && (s1.blue == s2.blue);
    }

    template<class T>
    bool operator!= (const RGBTmpl<T>& s1, const RGBTmpl<T>& s2)
    {
        return (s1.red != s2.red) || (s1.green != s2.green) || (s1.blue != s2.blue);
    }


    template<class T1, class T2>
    RGBTmpl<T1> operator+ (const RGBTmpl<T1>& s1, const T2& s2)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (s1.red   + s2);
        dst.green = Limit<T1> (s1.green + s2);
        dst.blue  = Limit<T1> (s1.blue  + s2);
        return dst;
    }
    // template <class S, class D> void Limit (const S& src, D& dst)

    template<class T1, class T2>
    RGBTmpl<T1> operator- (const RGBTmpl<T1>& s1, const T2& s2)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (s1.red   - s2);
        dst.green = Limit<T1> (s1.green - s2);
        dst.blue  = Limit<T1> (s1.blue  - s2);
        return dst;
    }

    template<class T1, class T2>
    RGBTmpl<T1> operator* (const RGBTmpl<T1>& s1, const T2& s2)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (s1.red   * s2);
        dst.green = Limit<T1> (s1.green * s2);
        dst.blue  = Limit<T1> (s1.blue  * s2);
        return dst;
    }

    template<class T1, class T2>
    RGBTmpl<T1> operator/ (const RGBTmpl<T1>& s1, const T2& s2)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (s1.red   / s2);
        dst.green = Limit<T1> (s1.green / s2);
        dst.blue  = Limit<T1> (s1.blue  / s2);
        return dst;
    }


    template<class T1>
    RGBTmpl<T1> operator+ (const RGBTmpl<T1>& s1, const double& d)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (d + s1.red);
        dst.green = Limit<T1> (d + s1.green);
        dst.blue  = Limit<T1> (d + s1.blue);
        return dst;
    }

    template<class T1>
    RGBTmpl<T1> operator- (const RGBTmpl<T1>& s1, const double& d)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (d - s1.red);
        dst.green = Limit<T1> (d - s1.green);
        dst.blue  = Limit<T1> (d - s1.blue);
        return dst;
    }

    template<class T1>
    RGBTmpl<T1> operator* (const RGBTmpl<T1>& s1, const double& d)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (d * s1.red);
        dst.green = Limit<T1> (d * s1.green);
        dst.blue  = Limit<T1> (d * s1.blue);
        return dst;
    }

    template<class T1>
    RGBTmpl<T1> operator/ (const RGBTmpl<T1>& s1, const double& d)
    {
        RGBTmpl<T1> dst;
        dst.red   = Limit<T1> (d / s1.red);
        dst.green = Limit<T1> (d / s1.green);
        dst.blue  = Limit<T1> (d / s1.blue);
        return dst;
    }


    template<class T>
    std::ostream& operator<< (std::ostream& stream, const RGBTmpl<T>& rgb)
    {
        stream << "(" << rgb.red << "," << rgb.green << "," << rgb.blue << ")";
        return stream;
    }

    // We want to display the pixel value as a number, not as a character
    // Even though this is an unsigned quantity, some compilers make mistakes
    template<> std::ostream& operator<< (std::ostream& stream, const RGBTmpl<Pel8>& rgb);



    //
    // RGBTmpl<T> mixed
    //
    //   Limit<RGBPel32s, RGB> (const RGBPel32s& src, RGB& dst);
    //
    //   3 parameter operations
    //   add2<D> (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   sub2<D> (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   mul2<D> (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   div2<D> (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //
    //   operator+ (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   operator- (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   operator* (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   operator/ (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //
    //   operator+= (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   operator-= (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   operator*= (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //   operator/= (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    //

    template<> RGB     Limit<RGB,     RGBPel32s> (const RGBPel32s& src);
    template<> RGB     Limit<RGB,     RGB>       (const RGB& src);
    template<> RGBPel32s Limit<RGBPel32s, RGBPel32s> (const RGBPel32s& src);


    template<class D, class S1, class S2>
    RGBTmpl<D> Min (const RGBTmpl<S1>& s1, const RGBTmpl<S2>& s2)
    {
        RGBTmpl<D> d;
        d.red   = Limit<D> (Min (s1.red,   s2.red));
        d.green = Limit<D> (Min (s1.green, s2.green));
        d.blue  = Limit<D> (Min (s1.blue,  s2.blue));
        return d;
    }


    template<class D, class S1, class S2>
    RGBTmpl<D> add2 (const RGBTmpl<S1>& s1, const RGBTmpl<S2>& s2)
    {
        RGBTmpl<D> d;
        d.red   = Limit<D> (s1.red   + s2.red);
        d.green = Limit<D> (s1.green + s2.green);
        d.blue  = Limit<D> (s1.blue  + s2.blue);
        return d;
    }

    template<class D, class S1, class S2>
    RGBTmpl<D> sub2 (const RGBTmpl<S1>& s1, const RGBTmpl<S2>& s2)
    {
        RGBTmpl<D> d;
        d.red   = Limit<D> (s1.red   - s2.red);
        d.green = Limit<D> (s1.green - s2.green);
        d.blue  = Limit<D> (s1.blue  - s2.blue);
        return d;
    }

    template<class D, class S1, class S2>
    RGBTmpl<D> mul2 (const RGBTmpl<S1>& s1, const RGBTmpl<S2>& s2)
    {
        RGBTmpl<D> d;
        d.red   = Limit<D> (s1.red   * s2.red);
        d.green = Limit<D> (s1.green * s2.green);
        d.blue  = Limit<D> (s1.blue  * s2.blue);
        return d;
    }

    template<class D, class S1, class S2>
    RGBTmpl<D> div2 (const RGBTmpl<S1>& s1, const RGBTmpl<S2>& s2)
    {
        RGBTmpl<D> d;
        d.red   = Limit<D> (s1.red   / s2.red);
        d.green = Limit<D> (s1.green / s2.green);
        d.blue  = Limit<D> (s1.blue  / s2.blue);
        return d;
    }

    template<class D, class S1>
    RGBTmpl<D> scale (const RGBTmpl<S1>& s1, float scaling)
    {
        RGBTmpl<D> d;
        d.red   = Limit<D> (scaling * static_cast<float>(s1.red));
        d.green = Limit<D> (scaling * static_cast<float>(s1.green));
        d.blue  = Limit<D> (scaling * static_cast<float>(s1.blue));
        return d;
    }


    template<class T1, class T2>
    RGBTmpl<T1> operator+ (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        RGBTmpl<T1> d;
        d.red   = Limit<T1> (s1.red   + s2.red);
        d.green = Limit<T1> (s1.green + s2.green);
        d.blue  = Limit<T1> (s1.blue  + s2.blue);
        return d;

        // On many platforms, a single line can be used, but some versions of
        // gcc have trouble figuring out which version of add2<> to use.
        // return add2<RGBTmpl<T1> > (s1, s2);
    }


    template<class T1, class T2>
    RGBTmpl<T1> operator- (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        RGBTmpl<T1> d;
        d.red   = Limit<T1> (s1.red   - s2.red);
        d.green = Limit<T1> (s1.green - s2.green);
        d.blue  = Limit<T1> (s1.blue  - s2.blue);
        return d;

        // See comment under operator+ (above)
        //return sub2<RGBTmpl<T1> > (s1, s2);
    }

    template<class T1, class T2>
    RGBTmpl<T1> operator* (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        RGBTmpl<T1> d;
        d.red   = Limit<T1> (s1.red   * s2.red);
        d.green = Limit<T1> (s1.green * s2.green);
        d.blue  = Limit<T1> (s1.blue  * s2.blue);
        return d;

        // See comment under operator+ (above)
        // return mul2<RGBTmpl<T1> > (s1, s2);
    }

    template<class T1, class T2>
    RGBTmpl<T1> operator/ (const RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        RGBTmpl<T1> d;
        d.red   = Limit<T1> (s1.red   / s2.red);
        d.green = Limit<T1> (s1.green / s2.green);
        d.blue  = Limit<T1> (s1.blue  / s2.blue);
        return d;

        // See comment under operator+ (above)
        // return div2<RGBTmpl<T1> > (s1, s2);
    }


    template<class T1, class T2>
    RGBTmpl<T1>& operator+= (RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        s1.red   = Limit<T1> (s1.red   + s2.red);
        s1.green = Limit<T1> (s1.green + s2.green);
        s1.blue  = Limit<T1> (s1.blue  + s2.blue);
        return s1;
    }

    template<class T1, class T2>
    RGBTmpl<T1>& operator-= (RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        s1.red   = Limit<T1> (s1.red   - s2.red);
        s1.green = Limit<T1> (s1.green - s2.green);
        s1.blue  = Limit<T1> (s1.blue  - s2.blue);
        return s1;
    }

    template<class T1, class T2>
    RGBTmpl<T1>& operator*= (RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        s1.red   = Limit<T1> (s1.red   * s2.red);
        s1.green = Limit<T1> (s1.green * s2.green);
        s1.blue  = Limit<T1> (s1.blue  * s2.blue);
        return s1;
    }

    template<class T1, class T2>
    RGBTmpl<T1>& operator/= (RGBTmpl<T1>& s1, const RGBTmpl<T2>& s2)
    {
        s1.red   = Limit<T1> (s1.red   / s2.red);
        s1.green = Limit<T1> (s1.green / s2.green);
        s1.blue  = Limit<T1> (s1.blue  / s2.blue);
        return s1;
    }

}

#endif // _imageTypes_h_

